% TODO: write the editor about this problem: "Color and shading may
% not be used" in tables,
% http://www.biomedcentral.com/bmcbioinformatics/authors/instructions/methodologyarticle#preparing-tables

% letter to the editor Dear editor, thank you for considering our
% article for publication. We have chosen to submit to BMC
% Bioinformatics because a couple of previous methods that we cite
% have appeared in BMC Bioinformatics, and we think the readers will
% appreciate the meaningful comparisons that we show between existing
% algorithms, using our new annotation-guided approach.

% I have spent much time trying to get the article to adhere to BMC
% Bioinformatics' "instructions for authors." I have just one concern
% with the editorial policy that "Color and shading may not be used"
% in tables, since I planned the article so that Figure 4 and Table 4
% appear on the same page, with a color code linking
% them. Additionally, the color code in Table 3 and Figure 3 is the
% same, and I feel that this greatly aides understanding and
% interpretation of our work. Is there any possibility that you can
% make an exception so that we can use color in the tables?



% Preparing a personal cover page (done, figure-smoothing-small.R and
% personal-cover-page-manual.png)

% If you wish to do so, you may submit an image which, in the event of
% publication, will be used to create a cover page for the PDF version
% of your article. The cover page will also display the journal logo,
% article title and citation details. The image may either be a figure
% from your manuscript or another relevant image. You must have
% permission from the copyright to reproduce the image. Images that do
% not meet our requirements will not be used.

% Images must be 300dpi and 155mm square (1831 x 1831 pixels for a
% raster image).

% Allowable formats - EPS, PDF (for line drawings), PNG, TIFF (for
% photographs and screen dumps), JPEG, BMP, DOC, PPT, CDX, TGF
% (ISIS/Draw).



% Reviewer's report
% Title: Learning smoothing models of copy number profiles using breakpoint
% annotations
% Version: 1 Date: 18 August 2012
% Reviewer number: 1
% Reviewer's report:

% Minor Essential Revisions

% The authors developed a tool for calling breakpoints in array-CGH
% data. They intended their tool to (i) accurately match the
% breakpoint calls made by an expert scientists by visual inspection,
% assuming that "the expert interpretation of the biologist is very
% accurate" (p2), and (ii) to be useable for data sets from different
% types of tumors after some initial training on a small number of
% sample profiles ("to generalize well to other data"; p2).

% That better support their claim that the newly developed tool indeed
% achieves these goals, the authors should

% 1. show additional, more complex DNA copy number profiles, including the
% breakponints called by the various tested models. The example shown in Fig. 1
% serves to illustrate the problem, but does not convincingly demonstrate the
% performance of their model. The breakpoints in Fig. 1 probably could be
% accurately called by any model as long as the parameters are somewhat
% reasonably chosen.

% 2. test their tool using at least one more data set of another tumor type.

% Level of interest: An article of importance in its field
% Quality of written English: Acceptable
% Statistical review: No, the manuscript does not need to be seen by a
% statistician.


% Reviewer's report
% Title: Learning smoothing models of copy number profiles using breakpoint
% annotations
% Version: 1 Date: 31 August 2012
% Reviewer number: 2
% Reviewer's report:
% Report on "Learning smoothing models of copy number profiles using
% breakpoint annotations" by Hocking et al.

% General

% The authors undertook a brave attempt to tackle an issue that I
% always thought will stay a subjective part of DNA copy number (CN)
% analysis: setting optimal smoothing parameters for breakpoint
% detection.  While I appraise this attempt, I do have some serious
% concerns which are addressed below.

% Major

% 1. Subjectivism is not completely eliminated, because subjective
% annotation of breakpoints is still required, although only for a
% limited number or regions and profiles. The annotator-bias may be
% serious [e.g.  I imagine that biologists may be tempted to find more
% breakpoints than a statistician, because they look differently at
% (noisy) data]. I would only trust this procedure if it would be
% based on multiple annotators, e.g. 3. The authors seem to realize
% this themselves, but chose not to address it ["In future work, it
% will be interesting to see if our conclusions are robust to the
% annotator"]. I think it is essential to address it. In fact, current
% practice of our microarray facility is to group-wise discuss profile
% plots from multiple breakpoint settings and come to a consensus on
% which one is best.

% 2. Is too many breakpoints as bad as too few? The authors seem to
% think so because they use a symmetric loss function. If the
% downstream analysis is based on breakpoint locations this may be the
% case. However, often the segmented data will be called (discretized)
% which may undo breakpoints of consecutive segments in the sense that
% the same discrete state is assigned to the segments.  Also, if the
% downstream analysis is based on the segment levels, introducing an
% extra segment (with level fairly close to its neighbor) does little
% harm. However, not detecting a breakpoint (as in Figure 1b,
% lambda=7.5, chromosome 7) may be more harmful.

% 3. What about small genomic regions? Are these as well observed by
% the annotators than large ones? Nowadays many biologists believe the
% small events to be more important than the large ones (in terms of
% driver and passenger alterations). This point is important when
% comparing algorithms. An algorithm like cghseg.k (with kmax=20)
% allows fairly few breakpoints, and hence the FPR is likely to be
% small. If the annotator does not well observe small genomic regions,
% also the "FNR" may be small; however FNR is underestimated due to
% the annotator's bias to observe breakpoints that are relatively far
% apart.

% 4. It seems the global and local model are two sub-optimal
% extremes. I'm sure that the authors are aware of shrinkage
% approaches, e.g. by pooling \lambda_i^p = p*\lambda_i + (1-p)
% \lambda, and finding optimal p in terms of mean square error. In
% particular in a large sample as the neuroblastoma one, this is
% likely to perform better than either of the two extremes. Why did
% they not pursue such an approach?

% Minor

% 1. p7. It seems min r_k and max r_k have not been defined.

% Level of interest: An article whose findings are important to those
% with closely related research interests

% Quality of written English: Acceptable
% Statistical review: Yes, and I have assessed the statistics in my report.


\NeedsTeXFormat{LaTeX2e}[1995/12/01]
\documentclass[10pt]{bmc_article}    


% TDH section titles
\newcommand{\sectionpick}{Selecting the optimal degree of smoothness}


\usepackage{cite}
\newcommand{\citep}[1]{\cite{#1}}
%\usepackage{natbib}
\usepackage{color} 
\usepackage{graphicx} 
\usepackage{amssymb,amsmath}

\usepackage{xcolor}
\usepackage{multirow}
\usepackage{textcomp}
%slashbox for tables with \ in upper left for row and column dimnames
%\usepackage{slashbox}
\usepackage{booktabs}
\newcommand{\url}[1]{\texttt{#1}}
\newcommand{\argmin}{\operatorname*{arg\, min}}
\newcommand{\model}[1]{#1}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\package}[1]{\texttt{#1}}
\input{algo.colors}

\newcommand{\NA}{\texttt{NA}}
\newcommand{\RR}{\mathbb R}
\newcommand{\NN}{\mathbb N}
% Load packages
\usepackage{url}  % Formatting web addresses  
\usepackage{ifthen}  % Conditional 
\usepackage{multicol}   %Columns
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails
\urlstyle{rm}
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%   
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %% 
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %% 
%%  submitted article.                         %% 
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                     


%\def\includegraphic{}
%\def\includegraphics{}



\setlength{\topmargin}{0.0cm}
\setlength{\textheight}{21.5cm}
\setlength{\oddsidemargin}{0cm} 
\setlength{\textwidth}{16.5cm}
\setlength{\columnsep}{0.6cm}

\newboolean{publ}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                              %%
%% You may change the following style settings  %%
%% Should you wish to format your article       %%
%% in a publication style for printing out and  %%
%% sharing with colleagues, but ensure that     %%
%% before submitting to BMC that the style is   %%
%% returned to the Review style setting.        %%
%%                                              %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

%Review style settings
\newenvironment{bmcformat}{\begin{raggedright}\baselineskip20pt\sloppy\setboolean{publ}{false}}{\end{raggedright}\baselineskip20pt\sloppy}

%Publication style settings
%\newenvironment{bmcformat}{\fussy\setboolean{publ}{true}}{\fussy}



% Begin ...
\begin{document}
\begin{bmcformat}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  \title{Learning smoothing models of copy number profiles using
    breakpoint annotations}
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Ensure \and is entered between all but   %%
%% the last two authors. This will be       %%
%% replaced by a comma in the final article %%
%%                                          %%
%% Ensure there are no trailing spaces at   %% 
%% the ends of the lines                    %%     	
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\author{Toby Dylan Hocking\correspondingauthor$^{1,2,3,4}$%
\email{Toby Dylan Hocking\correspondingauthor - toby.hocking@inria.fr}%
\and
Gudrun Schleiermacher$^5$%
\email{Gudrun Schleiermacher - gudrun.schleiermacher@curie.net}%
\and
Isabelle Janoueix-Lerosey$^{5}$%
\email{Isabelle Janoueix-Lerosey - Isabelle.Janoueix@curie.fr}%
\and
Valentina Boeva$^{4}$%
\email{Valentina Boeva - Valentina.Boeva@curie.fr}%
\and
Julie Cappo$^{5}$%
\email{Julie Cappo - Julie.Cappo@curie.fr}%
\and
Olivier Delattre$^5$%
\email{Olivier Delattre - Olivier.Delattre@curie.fr}%
\and
Francis Bach$^1$%
\email{Francis Bach - francis.bach@ens.fr}
and
Jean-Philippe Vert$^{2,3,4}$%
\email{Jean-Philippe Vert - Jean-Philippe.Vert@mines-paristech.fr}%
}
      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address{%
  \iid(1) INRIA -- Sierra project-team, Paris, F-75013, France\\
  \iid(2) Centre for computational biology, Mines ParisTech,
  Fontainebleau, F-77300, France\\
  \iid(3) Institut Curie, Paris, France\\
  \iid(4) INSERM U900, Paris, F-75248, France\\
  \iid(5) INSERM U830, Paris, F-75248, France}%

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% The Section headings here are those for  %%
%% a Research article submitted to a        %%
%% BMC-Series journal.                      %%  
%%                                          %%
%% If your article is not of this type,     %%
%% then refer to the Instructions for       %%
%% authors on http://www.biomedcentral.com  %%
%% and change the section headings          %%
%% accordingly.                             %%   
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{abstract}
        % Do not use inserted blank lines (ie \\) until main body of text.
        \paragraph*{Background:} 
Many models have been proposed to detect copy number alterations in
chromosomal copy number profiles, but it is usually not obvious to
decide which is most effective for a given data set. Furthermore, most
methods have a smoothing parameter that determines the number of
breakpoints and must be chosen using various heuristics.

\paragraph*{Results:} 
We present three contributions for copy number profile smoothing model
selection. First, we propose to select the model and degree of
smoothness that maximizes agreement with visual breakpoint region
annotations. Second, we develop cross-validation procedures to
estimate the error of the trained models. Third, we apply these
methods to a new database of annotated neuroblastoma copy number
profiles, which we make available as a public benchmark for testing
new algorithms.

\paragraph*{Conclusions:} 
We created 2 databases of breakpoint annotations for a set of 575
neuroblastoma copy number profiles, and tested 17 smoothing models to
determine the most accurate breakpoint detector in these data.
Whereas previous studies have been qualitative or limited to simulated
data, our approach is quantitative and suggests which algorithms are
fastest and most accurate in practice on real data.

\end{abstract}



\ifthenelse{\boolean{publ}}{\begin{multicols}{2}}{}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% The Section headings here are those for  %%
%% a Research article submitted to a        %%
%% BMC-Series journal.                      %%  
%%                                          %%
%% If your article is not of this type,     %%
%% then refer to the instructions for       %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and change the section headings          %%
%% accordingly.                             %% 
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%
%% Background %%
%%
\section*{Background: the need for smoothing model selection
  criteria}

DNA copy number alterations (CNAs) can result from various types of
genomic rearrangements, and are important in the study of many types
of cancer \citep{weinberg}. In particular, clinical outcome of
patients with neuroblastoma has been shown to be worse for tumors with
segmental alterations or breakpoints in specific genomic regions
\citep{isabelle-2009,gudrun-jclinicaloncology}.  Thus, to construct an
accurate predictive model of clinical outcome for these tumors, we
must first accurately detect the precise location of each breakpoint.


In recent years, array comparative genomic hybridization (aCGH)
microarrays have been developed as genome-wide assays for CNAs, using
the fact that microarray fluoresence intensity is proportional to DNA
copy number \citep{pinkel}. In parallel, there have been many new
mathematical models proposed to smooth the noisy signals from these
microarray assays in order to recover the CNAs
\citep{glad,statistical-approach,dnacopy,cghFLasso,haarseg,fused-lasso-path,gada,pelt}.
Each model has different assumptions about the data, and it is not
obvious to decide which model is appropriate for a given data set.

Furthermore, most models have parameters that control the degree of
smoothness. Varying these smoothing parameters will vary the number of
detected breakpoints. Most authors give default values that accurately
detect breakpoints on some data, but do not necessarily generalize
well to other data. There are some specific criteria for choosing the
degree of smoothness in some models
\citep{lavielle2005,mBIC,penalized-cna}, but it is impossible to
verify whether or not the mathematical assumptions of these models are
satisfied for real noisy microarray data.

To motivate the use of their cghFLasso smoothing model,
Tibshirani and Wang write ``The results of a CGH experiment are often
interpreted by a biologist, but this is time consuming and not
necessarily very accurate'' \cite{cghFLasso}.
% quote from penncnv article, Parameter estimation and CNV calling
% "The initial model parameters for the HMM were estimated empirically
% from several large CNV regions, through manually examining the
% BeadStudio Genome Viewer for a set of genotyped individuals."
In contrast, this paper takes the opposite view and assumes that the
expert interpretation of the biologist is the gold standard which a
model should attain. The first contribution of this paper is a
smoothing model training protocol based on this assumption.

In practice, visualization tools such as VAMP \citep{vamp} are used to
plot the normalized microarray signals against genomic position for
interpretation by an expert biologist looking for CNAs. Then the
biologist plots a model and varies its smoothness parameter, until the
model seems to capture all the visible breakpoints the data. In this
article, we make this model training protocol concrete by using an
annotation database to encode the expert's interpretation.

The particular type of annotations that we propose are counts of
breakpoints in genomic regions. By visual inspection of the noisy
signal, it is not obvious to locate the exact location of a
breakpoint, but it is easy to determine whether or not a region
contains a breakpoint. So rather than defining annotations in terms of
precise breakpoint locations, we instead define them in terms of
regions. For every region, we record the number of breakpoints that an
expert expects in that region. These annotated regions can then be
used to select an appropriate model, as shown in Figures~1 and 2.

We note that using databases of visual annotations is not a new idea,
and has been used successfully for object recognition in photos and
cell phenotype recognition in microscopy
\citep{labelme,cellprofiler}. In array CGH analysis, some models can
incorporate prior knowledge of locations of CNAs \cite{shah}, but no
models have been specifically designed to exploit visual breakpoint
annotations.

Our second contribution is a protocol to estimate the breakpoint
detection ability of the trained smoothing models on real data. In the
Methods section, we propose to estimate the false positive and false
negative rates of the trained models using cross-validation. This
provides a quantitative criterion for deciding which smoothing
algorithms are appropriate breakpoint detectors for which data.

The third contribution of this paper is a systematic, quantitative
comparison of the accuracy of 17 common smoothing algorithms on a
new database of 575 annotated neuroblastoma copy number profiles,
which we give in the Results section. There are several publications
which attempt to assess the accuracy of smoothing algorithms, and
these methods fall into 2 categories: simulations and low-throughput
experiments. GLAD, DNAcopy, and a hidden Markov model were compared by
examining false positive and false negative rates for detection of a
breakpoint at a known location in simulated data
\cite{compare}. However, there is no way to verify if the assumptions
of the simulation hold in a real data set, so the value of the
comparison is limited. In another article, the accuracy of the
\model{CNVfinder} algorithm was assessed using quantitative PCR
\cite{cnvfinder}. But quantitative PCR is low-throughput and costly,
so is not routinely done as a quality control. So in fact there are no
previous studies that quantitatively compare breakpoint detection of
smoothing models on real data. In this paper we propose to use
annotated regions instead of precise breakpoint locations for
quantifying smoothing model accuracy, and we make available 575 new
annotated neuroblastoma copy number profiles as a benchmark for the
community to test new algorithms on real data.

Several authors have recently proposed methods for so-called joint
segmentation of multiple CGH profiles, under the hypothesis that each
profile shares breakpoints in the exact same location
\citep{jp-nips,nbc}. These models are not useful in our
setting, since we assume that breakpoints do not occur in the exact
same locations across copy number profiles. Instead, we focus on the
case of learning a model that will accurately detect an unknown number
of breakpoints in a copy number profile.

To summarize, this article describes a quantitative method for DNA
copy number profile smoothing model selection. First, an expert
examines scatterplots of the data, and encodes her interpretation of
the breakpoint locations in a database of annotated regions. To
repeat, the annotations represent an expert's interpretation, not the
biological truth in the tumor samples, which is unknown. We treat the
annotated regions as a gold standard, and compare them to the
breakpoints detected by 17 existing models. The best model for our
expert is the one which maximizes agreement with the annotation
database.


\section*{Results and discussion}\label{results}

The 17 smoothing models described in the Algorithms section were
applied to 575 neuroblastoma copy number profiles, described in the
Data section. Note that to decrease computation time, the model
fitting may be trivially parallelized for profiles, algorithms, and
smoothing parameter values.

After fitting the models, we used breakpoint annotations to quantify
the accuracy of each model. We constructed 2 annotation databases
based on 2 different experts' interpretations of the same 575 profiles
(Table~1). The ``original'' annotations were created by typing 0 or 1
in a 6-column spreadsheet after systematic inspection of the same 6
regions on each profile. The ``detailed'' annotations were constructed
by using GUIs which allow zooming and direct annotation on the plotted
profiles. The 2 annotation data sets are mostly consistent, but the
detailed annotations provide more precise breakpoint locations
(Figure~3).

For both annotation databases, we calculated the global and local
error curves $E^{\text{global}}(\lambda)$ and
$E_i^{\text{local}}(\lambda)$, which quantify how many breakpoint
annotations disagree with the model breakpoints. As shown in Figure~2
for original set of annotations, the smoothness parameter $\lambda$ is
chosen by minimizing the error curves.


\subsection*{Among global models, \model{cghseg.k} 
and \model{pelt.n}
  exhibit the smallest training error}

The global model is defined as the smoothness parameter $\hat \lambda$
that minimizes the global error $E^{\text{global}}(\lambda)$, which is
the total number of incorrect annotations over all profiles.  Training
error curves for cghseg.k, pelt.n, flsa.norm, dnacopy.sd, and
glad.lambdabreak are shown in Figure~2. An ideal global model would
have zero annotation error $E^{\text{global}}(\hat \lambda)=0$ for
some smoothness parameter $\hat \lambda$. However, none of the 17
global models that we examined achieved zero training error in either
of the two annotation databases. The best global models were the
equivalent \model{cghseg.k} and \model{pelt.n} models, which achieved
the minimum error of \input{global-original}\unskip\% and
\input{global-detailed}\unskip\% in the original and detailed data
sets.


The ROC curves for the training error of the global models for each
algorithm are traced in Figure~4. It is clear that the default
parameters of each algorithm show relatively large false positive
rates. The only exception is the pelt.default algorithm, which showed
low false positive and true postive rates. The models chosen by
maximizing agreement with the breakpoint annotation data also exhibit
smaller false positive rates at the cost of smaller true positive
rates. The ROC curves suggest that the equivalent \model{cghseg.k} and
\model{pelt.n} models are the most discriminative for breakpoint
detection in the neuroblastoma data.

\subsection*{Among local models, \model{cghseg.k} and \model{pelt.n}
  exhibit the smallest training error}

Since there is no global model that agrees with all of the annotations
in either database, we fit local models with profile-specific
smoothness parameters $\hat \lambda_i$. For every profile $i$, the
local model is defined as the smoothness parameter $\hat \lambda_i$
that minimizes the local error $E_i^{\text{local}}(\lambda)$, the
number of incorrect annotations with respect to the annotations on
profile $i$. As shown in Figure~2, the local model fits the
annotations at least as well as the global model:
$E_i^{\text{local}}(\hat \lambda_i)\leq E^{\text{global}}(\hat
\lambda)$. However, the local model does not necessarily attain zero
error. For example, Figure~2 shows that \model{dnacopy.sd} does not
detect a breakpoint in profile $i=362$ even at the smallest parameter
value, corresponding to the model with the most breakpoints.

In Table~2, we compare the fitting ability of the local models on the
2 annotation data sets. It clearly shows that some models are better
than others for fitting the expert annotations of the neuroblastoma
data. In particular, the equivalent \model{cghseg.k} and
\model{pelt.n} local models show the best fit, with
\input{local-original}\unskip\% and \input{local-detailed}\unskip\%
error in the original and detailed data sets. Note that these are
lower error rates than the global models, as expected.

But even if the local models are better at fitting the given
breakpoint annotations, they do not generalize well to un-annotated
breakpoints, as we show in the next section.

\subsection*{Global models detect un-annotated breakpoints 
better than
  default models}

Leave-one-out cross-validation was used to estimate the breakpoint
detection of each model. Figure~5 shows the error rates of each model,
across both annotation data sets.  It is clear that the training
procedure makes no difference for models pelt.default, glad.default,
dnacopy.default, cghseg.mBIC, gada.default and cghFLasso, which are
default models with no smoothness parameters. Each of these models is
inferior to its respective global model in terms of breakpoint
detection. The large error of these models suggest that the
assumptions of their default parameter values do not hold in the
neuroblastoma data set. More generally, these error rates suggest that
smoothness parameter tuning is critically important to obtain an
accurate smoothing of real copy number profiles.

To show an example of how the learned models outperform default
models, Figure~6 shows one representative profile with many
breakpoints. Note that the models were trained on other profiles, so
the shown annotations can be used for model evaluation. For this
profile, dnacopy.default shows 2 false positives and 2 false
negatives, and dnacopy.sd shows no improvement with 4 false
negatives. The pelt.default and cghseg.mBIC show 10 and 3 annotation
errors, respectively. The cghseg.k and pelt.n global models show only
2 annotation errors, demonstrating the usefulness of annotation-based
model training.

In addition, there are 5 supplementary figures which compare these
default and global models. Again, the global models were learned on
other profiles, so the shown annotations can be used for model
evaluation.

\subsection*{Global models detect un-annotated
 breakpoints better than local
  models}

The leave-one-out cross-validation results in Figure~5 also allow
comparison of global and local models. For dnacopy.prune,
glad.MinBkpWeight, glad.lambdabreak, dnacopy.sd and flsa, there is
little difference between the local and global training
procedures. For models flsa.norm, gada, pelt.n, and cghseg.k, there
is a clear advantage for the global models which share
information between profiles. The equivalent cghseg.k and pelt.n
models show the minimal test error of only
\input{global-test-original}\unskip\% and
\input{global-test-detailed}\unskip\% in the original and detailed
data sets.

% Finally, note that the false positive rate of locally trained models
% is higher than the false negative rate for most algorithms. This can
% be explained by the larger fraction of normal annotations present in
% the training set, and the fact that many profiles have only normal
% annotations (Table~2).

\subsection*{Only a few profiles need to be annotated for a good global
  model}

To estimate the generalization error of a global model trained on a
relatively small training set of $t$ annotated profiles, we applied
$\lfloor n/t\rfloor $-fold cross-validation to the $n=575$ profiles.

For several training set sizes $t$, we plot the accuracy of the
cghseg.k, pelt.n, gada, flsa.norm, dnacopy.sd and glad.lambdabreak
models in Figure~7. It shows that adding more annotations to the
training set increases the breakpoint detection accuracy in general,
but at a diminishing rate. Each model quickly attains its specific
maximum, after only about $t=10$ training profiles.

In Table~3, we used $\lfloor n/t\rfloor $-fold cross-validation in the
detailed annotations to estimate the error rates of all 17 models
trained using only $t=10$ profiles.
%These error estimates are slightly larger, but the model
%ordering is mostly unchanged, with respect to the leave-one-out
%cross-validated estimates of the global model error rates in Table~3.
The equivalent cghseg.k and pelt.n models show the best performance on
these data, with an estimated breakpoint detection error of $ 7.7\%$.

\subsection*{Global models generalize across annotators}

We assessed the extent to which the annotator affects the results by
comparing models trained on one data set and tested on the
other. Figure~7 shows that test error changes very little between
models trained on one data set or the other. This demonstrates that
global models generalize very well across annotators.

\subsection*{Timing PELT and cghseg}
The PELT and cghseg models use different algorithms to calculate the
same segmentation, which showed the best breakpoint detection
performance in every comparison. But they are slightly different in
terms of speed, as we show in Table~3.

When comparing the global models, cghseg.k is somewhat faster than
pelt.n. For cghseg.k, pruned dynamic programming is used to calculate
the best segmentation $\mu^k$ for $k\in\{1,\dots,20\}$ segments, which
is the slow step. Then, we calculate the best segmentation for
$\lambda\in\{\lambda_1,\dots,\lambda_{100}\}$, based on the stored
$\mu^k$ values. In contrast, the Pruned Exact Linear Time algorithm
must be run for each $\lambda\in\{\lambda_1,\dots,\lambda_{100}\}$,
and there is no information shared between $\lambda$ values.

Timing the PELT and cghseg default models without tuning
parameters shows the opposite trend. In particular, the default
cghseg.mBIC method is slower than the pelt.default method. This makes
sense since cghseg must first calculate the best segmentation $\mu^k$
for several $k$, then use the mBIC criterion to choose among them. In
contrast, the PELT algorithm recovers just the $\mu^k$ which
corresponds to the Schwarz Information Criterion penalty constant
$\beta=\log d$. So if you want to use a particular penalty constant
$\beta$ instead of the annotation-guided approach we suggest in this
chapter, the default PELT method offers a modest speedup over cghseg.

\subsection*{Annotation-based modeling is feasible for high-density
  microarrays}

Although not the main focus of this paper, we have already started to
apply annotation-based modeling to high-density microarrays. For
example, Figure~8 shows part of chromosome 2 from an Affymetrix SNP6
microarray. This microarray offers almost 2 million probes, and we
show annotated breakpoints around 3 CNAs from $\approx 1$Mb to
$\approx 10$kb. As long as there is GUI software that supports zooming
and annotation, it is feasible to apply annotation-based modeling.


\section*{Conclusions}

We proposed to train breakpoint detection models using annotations
determined by visual inspection of the copy number profiles. We have
demonstrated that this approach allows quantitative comparison of
smoothing models on a new data set of 575 neuroblastoma copy number
profiles. These data provide the first set of annotations that can be
used for benchmarking the breakpoint detection ability of future
algorithms. Our annotation-based approach is quite useful in practice
on real data, since it provides a quantitative criterion for choosing
the model and its smoothing parameter. 
%We showed that learning a
%global smoothness parameter on a limited set of annotations can
%generalize well to un-annotated profiles. 

One possible criticism of annotation-based model selection is the time
required to create the annotations. However, using the GUIs that we
have developed, it takes only a few minutes to annotate the
breakpoints in a profile. This is a relatively small investment
compared to the time required to write the code for data analysis,
which is typically on the order of days or weeks. In addition, in the
neuroblastoma data, we observed that annotating only about 10 of 575
profiles was sufficient to learn a smoothness parameter that achieves
the model-specific optimal breakpoint detection. More generally, our
results suggest that after obtaining a moderately sized database of
annotations, data analysis time is better spent designing and testing
better models. Additionally, the learned models generalized very well
between annotators. So breakpoint annotations are a feasible approach
for finding an accurate model and smoothing parameter for real copy
number profiles.


% Furthermore, we propose to streamline the annotation process by
% using a simple, portable, free software Python annotation GUI which is
% available as package \texttt{annotate\_regions} from the Python
% Package Index.

% Another potential criticism of our approach is that annotations
% determined by visual inspection may not be consistent between
% experts. To investigate this concern in future work, we could build a
% database of several experts' annotations of the same profiles. In any
% case, it is clear that our approach will favor algorithms that capture
% breakpoints that agree with the annotator's visual definition of a
% breakpoint.

We compared local models for single profiles with global models
selected using annotations from several profiles. We observed that
local models fit the given annotations better, but global models
generalize better to un-annotated regions. In contrast with our
results, it has been claimed that local models should be better in
some sense: ``it is clear that the advantages of selecting
individual-specific $\lambda$ values outweigh the benefit of selecting
constant $\lambda$ values that maximize overall performance''
\cite{penalized-cna}.  However, they did not demonstrate this claim
explicitly, and one of the contributions of this work is to show that
global models generalize better than local models, according to our
leave-one-out estimates. 


% However, the smoothness parameterization must be
% carefully chosen. For example, the flsa.norm algorithm scales the
% smoothness parameter $\lambda$ by the number of points and the length
% of the chromosome, and results in lower error rates than the unscaled
% flsa algorithm. We are interested in investigating other
% parameterizations which could reduce the error even further.



It will be interesting to apply annotation-based model training to
other algorithms and data sets. In both annotation data sets we
analyzed, \model{cghseg.k} and \model{pelt.n} showed the best
breakpoint detection, but another model may be selected for other
data.
 
Our results indicate that even the best models have non-zero training
and testing breakpoint detection error, which could be
improved. To make a model that perfectly fits the training
annotations, a dynamic programming algorithm called SegAnnot was
proposed to recover the most likely breakpoints that are consistent
with the annotation data \cite{SegAnnot}. Developing a model that
lowers the test error remains an interesting direction of future
research.


We have solved the problem of smoothness parameter selection using
breakpoint annotations, but the question of detecting CNAs remains. By
constructing a database of annotated regions of CNAs, we could use a
similar approach to train models that detect CNAs. Annotations could
be actual copy number ($0,1,2,3,\dots$) or some simplification (loss,
normal, gain).  We will be interested in developing joint breakpoint
detection and copy number calling models that directly use these
annotation data as constraints or as part of the model likelihood.


\section*{Methods}
\label{training}
% The first contribution of this paper is a smoothing model training
% protocol that uses visually determined breakpoint annotations to
% quantify model accuracy, which we explain in this section.

\subsection*{GUIs for annotating copy number profiles}

Assume that we have $n$ DNA copy number profiles, and we would like to
accurately detect their breakpoints. The first step of
annotation-based modeling is to plot the data, visually identify
breakpoints, and save these regions to an annotation database. We
created 2 annotation GUIs for this purpose: a Python program
for low-density profiles called
\mbox{\texttt{annotate\_breakpoints.py}}, and a web site for larger
profiles called SegAnnDB.

We used Tkinter in Python's standard library to write
\texttt{annotate\_breakpoints.py}, a cross-platform GUI for annotating
low-density DNA copy number profiles. The annotator loads several
profiles from a CSV file, plots the data, and allows annotated regions
to be drawn on the plot and saved to a CSV file for later
analysis. The annotator does not support zooming so is not suitable
for annotating high-density profiles. It is available in the
\texttt{annotate\_regions} package on the Python Package Index:

\url{http://pypi.python.org/pypi/annotate_regions}

SegAnnDB is a web site that can be used to annotate low to
high-density copy number profiles. After copy number data in bedGraph
format is uploaded, the server uses D3 to show plots which can be
annotated \citep{d3}. The annotations can then be downloaded for later
analysis. As shown in Figure~8, the plots can be zoomed for detailed
annotation of high-density copy number profiles. The free/open-source
software that runs the web site can be downloaded from the
\texttt{breakpoints} project on INRIA GForge:

\url{https://gforge.inria.fr/scm/viewvc.php/webapp/modeler/?root=breakpoints}

% And a prototype of the annotation web site is online:

% \url{https://bioviz.rocq.inria.fr/plotter/}

\subsection*{Definition of breakpoints in smoothing models}

%Assume we have $n$ annotated DNA copy number profiles. 
For each profile $i\in\{1,\dots,n\}$, we observe $d_i\in\NN$ noisy
logratio measurements $y_i\in\RR^{d_i}$. Assume that we have a model
with smoothness parameter $\lambda$ that takes the vector of logratios
$y_i$ and outputs a smoothed signal $\hat
y_i^\lambda\in\RR^{d_i}$. For simplicity of notation, let $\hat
x^\lambda\in\RR^m$ be the smoothed signal sampled at positions $p_1\leq \dots
\leq p_m$ on one chromosome of one profile. We define the breakpoints
predicted by this model as the set of positions where there are jumps
in the smoothed signal:
\begin{equation}
  \label{eq:predicted}
  \hat b^\lambda = \big\{
(p_j+p_{j+1})/2
\mid
\hat x_{j}^\lambda\neq
\hat x_{j+1}^\lambda,
\ \forall j\in\{1,\dots,m-1\}
\big\}
\end{equation}
Note that this set is drawn using vertical black lines in Figures~1
and 6.

\subsection*{Definition of the annotation error}

For every profile and chromosome, we judge the accuracy of the
predicted set of breakpoints $\hat b^\lambda$ using a set of
visually-determined regions and corresponding annotations. Every
annotation $a=[\underline a,\overline a]$ is an interval that
specifies the expected number of changes in the corresponding region
$r$. For example, we defined 3 types of annotations: $a=[0,0]$ for
0breakpoints annotations, $a=[1,1]$ for 1breakpoint annotations, and
$a=[1,\infty)$ for $>$0breakpoints annotations. A region is an
interval of base pairs on the chromosome that corresponds to $\hat
b^\lambda$, for example $r=[1000000,2000000]$.

The false positives (FP) and false negatives (FN) are calculated by
comparing the estimated number of changes in each region $|\hat
b^\lambda\cap r|$ to the annotated number of changes $a$ using the
zero-one loss:
\begin{equation}
  \label{eq:FP}
 \text{FP}(\hat b^\lambda, r, a)
 = 
%\sum_{(r,a)\in(R,A)}
 \begin{cases}
   1 & \text{ if }|\hat b^\lambda\cap r|> \overline a,\\
   0 & \text{ otherwise.}
 \end{cases}
\end{equation}
\begin{equation}
  \label{eq:FN}
 \text{FN}(\hat b^\lambda, r, a)
 = 
%\sum_{(r,a)\in(R,A)}
 \begin{cases}
   1 & \text{ if }|\hat b^\lambda\cap r|< \underline a,\\
   0 & \text{ otherwise.}
 \end{cases}
\end{equation}
Note that $>$0breakpoints annotations can never have false positives
$\text{FP}(\lambda, r, [1,\infty) )=0$ and 0breakpoints annotations
can never have false negatives $\text{FN}(\lambda, r, [0,0])=0$. The
annotation error is defined as the sum of false positives and false
negatives:
\begin{equation}
  \label{eq:e}
  e(\hat b^\lambda, r, a) = \text{FP}(\hat b^\lambda, r, a)+
  \text{FN}(\hat b^\lambda, r, a)=
 \begin{cases}
   0 & \text{ if }|\hat b^\lambda\cap r|\in a\\
   1 & \text{ otherwise.}
 \end{cases}
\end{equation}
Note that this loss function gives the same weight to false positives
and false negatives. Re-weighting schemes could be used, but uniform
weighting is justified in the data we analyzed since each annotation
took approximately the same amount of time to create.

\subsection*{Definitions of error and ROC curves}

For each profile $i$, we define the local error as the total
annotation error over all annotated regions:
\begin{equation}
  \label{eq:localerror}
  E_i^{\text{local}}(\lambda)=\sum_{(\hat b^\lambda, r, a)\text{ on profile $i$}} 
e(\hat b^\lambda, r, a).
\end{equation}
We define the global error as the total annotation error over all
profiles:
\begin{equation}
  \label{eq:globalerror}
  E^{\text{global}}(\lambda)=\sum_{i=1}^n E_i^{\text{local}}(\lambda).
\end{equation}

For a given algorithm, a ROC curve is drawn by plotting the true
positive rate $\text{TPR}(\lambda)$ against the false positive rate
$\text{FPR}(\lambda)$ for all values of $\lambda$. For one
annotation, the true positive indicator function is
\begin{equation}
  \label{eq:tp}
  \text{TP}(\hat b^\lambda, r, a) =
  \begin{cases}
   1 & \text{ if }|\hat b^\lambda\cap r|\geq \underline a\\
   0 & \text{ otherwise.}
  \end{cases}
\end{equation}
To define the true positive rate, we first define the set of positive
annotations $(\hat B^\lambda_+,R_+,A_+)$ as all the annotations $a\in
A_+$ such that there is at least one breakpoint $\underline a\geq
1$. The true positive rate over all positive annotations is
\begin{equation}
  \label{eq:tpr}
  \text{TPR}(\lambda) = 
\frac{1}{|A_+|}
\sum_{(\hat b^\lambda, r, a)\in (\hat B^\lambda_+,R_+,A_+)}
\text{TP}(\hat b^\lambda, r, a).
\end{equation}
To define the false positive rate, we first define the set of negative
annotations $(\hat B^\lambda_-, R_-, A_-)$ as all the annotations that
could possibly have a false positive: $\overline a<\infty$. Then the
false positive rate over all negative annotations is
\begin{equation}
  \label{eq:tpr}
  \text{FPR}(\lambda) = 
\frac{1}{|A_-|}
\sum_{(\hat b^\lambda, r, a)\in (\hat B^\lambda_-,R_-,A_-)}
\text{FP}(\hat b^\lambda, r, a).
\end{equation}

\label{quantify}

\subsection*{\sectionpick}

We assume that $\lambda$ is a tuning parameter that is monotonic in
the number of breakpoints, which is the case for the models considered
in this paper. Fix a set of smoothing parameters $\lambda\in\Lambda$,
and run the smoothing algorithm with each of these
parameters. Intuitively, we should select the value of $\lambda$ that
maximizes agreement with annotation data. For global models, we
minimize the global error, and there is usually one best value:
\begin{equation}
  \label{eq:hat_lambda_i}
  \hat \lambda=\argmin_{\lambda\in\Lambda}\  E^\text{global}(\lambda).
\end{equation}

For the local model for profile $i$, we want to minimize the local
error:
\begin{equation}
  \label{eq:hat_lambda_i}
  \hat \lambda_i=\argmin_{\lambda\in\Lambda}\  E_i^\text{local}(\lambda).
\end{equation}
Since the training set consists of only the annotations of one profile
$i$, there may be several smoothing parameters $\lambda$ that minimize
the error.
\label{pick} We propose to choose between models that achieve the
minimum error based on the shape of the error curve, and these cases
are illustrated in Figure~2.
% \begin{equation}
%   \label{eq:pick}
%   L(E)=
%   \begin{cases}
%     \lambda^* & \text{ if }E(\lambda_1)=E(\lambda_\ell)=E^*\\
%     \lambda^*_{\ell^*}&\text{ if }E(\lambda_1)=E^*\\
%     \lambda^*_{1} &\text{ if }E(\lambda_\ell)=E^*\\
%     \lambda^*_{\lceil \ell^*/2\rceil} & \text{ otherwise}
%   \end{cases}
% \end{equation}
\begin{enumerate}
\item When the minimum error is achieved in a range of intermediate
  parameter values, we select a value in the middle. This occurs in the
  local error curves shown for \model{flsa.norm} and
  \model{cghseg.k}.
\item When the minimum is attained by the model with the most
  breakpoints, we select the model with the fewest breakpoints that has
  the same error. This attempts to minimize the false positive
  rate. This occurs for profile $i=375$ with
  the \model{dnacopy.sd} model.
\item When the minimum is attained by the model with the fewest
  breakpoints, we select the model with the most breakpoints that has
  the same error. This attempts to minimize the false negative rate,
  and occurs for profile $i=362$ with 
  the \model{dnacopy.sd} model.
\end{enumerate}
More complicated smoothing parameter estimators could be defined, but
for simplicity in this article we explore only the global $\hat
\lambda$ and local $\hat \lambda_i$ models.


\subsection*{Leave-one-out cross-validation for comparing local and
  global models}
To compare the breakpoint detection performance of local and global
models, we propose to leave one annotation per profile aside as a test
set. The input parameter $V$ is the number of times the procedure is
repeated. In our analysis we took $V=10$ repetitions. For each
repetition,
\begin{enumerate}
\item On each profile, randomly pick one annotated region and set it
  aside as a test set.
\item Using all the other annotations as a training set, select the
  best $\lambda$ using the protocol described in
%Section~\ref{pick}. 
  Section~``\sectionpick.''  For local models we learn a
  profile-specific $\hat \lambda_i$ that minimizes
  $E_i^{\text{local}}$, and for global models we learn a global $\hat
  \lambda$ that minimizes $E^{\text{global}}$.
\item To estimate how the model generalizes, count the errors of the
  learned model on the test regions.
\end{enumerate}
The final estimate of model error shown in Figure~5 is the average
error over all $V$ repetitions.

\subsection*{$\lfloor n/t\rfloor$-fold cross-validation to estimate
  error on un-annotated profiles}

Since the annotation process is time-consuming, we are interested in
training an accurate breakpoint detector with as few annotations as
possible. Thus we would like to answer the following question: how
many profiles $t$ do I need to annotate before I get a global model
that will generalize well to all the other profiles?

To answer this question, we estimate the error of a global model
trained on the annotations from $t$ profiles using
cross-validation. We divide the set of $n$ annotated profiles into
exactly $\lfloor n/t \rfloor$ folds, each with approximately $t$
profiles. For each fold, we consider its annotations a training set
for a global model, and combine the other folds as a test set to
quantify the model error. The final estimate of generalization error
is then the average model error over all folds.



\subsection*{Data: neuroblastoma copy number profiles}\label{data}

We analyzed a new data set of $n=575$ copy number profiles from aCGH
microarray experiments on neuroblastoma tumors taken from patients at
diagnosis. The microarrays were produced using various technologies,
so do not all have the same probes. The number of probes per
microarray varies from 1719 to 71340, and the median spacing between
probes varies from 40Kb to 1.2Mb. In this article we analyzed the
normalized logratio measurements of these microarrays, which we have
made available as \code{neuroblastoma\$profiles} in R package
\package{neuroblastoma} on CRAN.

Two different expert annotations were used to construct 2 annotation
databases based on these profiles (Table~1). The ``original''
annotations were created by typing 0 or 1 in a 6-column spreadsheet
after systematic inspection of the same 6 regions on each
profile. These annotations are available as
\code{neuroblastoma\$annotations} in R package \package{neuroblastoma}
on CRAN. The ``detailed'' annotations were constructed by using GUIs
which allow zooming and direct annotation on the plotted profiles. The
2 annotation data sets are mostly consistent, but the detailed
annotations provide more precise breakpoint locations (Figure~3).
These annotations are available as \code{neuroblastomaDetailed} in R
package \package{bams} on CRAN.

\subsection*{Algorithms: copy number profile smoothing
  models}\label{models}

In this study we considered smoothing models from the bioinformatics
literature with free software implementations available as R packages
on CRAN, R-Forge, or Bioconductor \citep{R,R-Forge,Bioconductor}. For
each algorithm, we considered three types of training for the
smoothness parameter $\lambda$:
\begin{itemize}
\item \textbf{Default models} can be used when functions give default
  parameter values, or do not have smoothness parameters that vary the
  number of breakpoints.
\item \textbf{Local models} choose a smoothness parameter that
  maximizes agreement with the annotations from a single profile.
\item \textbf{Global models} choose a smoothness parameter
  that maximizes agreement with the entire database of training
  annotations.
\end{itemize}
In the following paragraphs, we discuss the precise meaning of the
smoothness parameter $\lambda$ in each of the algorithms. The code
that standardizes the outputs of these models can be found in the list
of functions \code{smoothers} in R package \package{bams} on CRAN. For
some algorithms (GADA, GLAD, DNAcopy) the smoothing $\hat y^\lambda\in\RR^d$
is defined for an entire profile $y\in\RR^d$, but in others (cghseg,
pelt, flsa) the smoothing $\hat x^\lambda\in\RR^m$ is defined in terms
of probes on a single chromosome $x\in\RR^m$.

We used version 1.0 of the \package{gada} package from R-Forge to
calculate a sparse Bayesian learning model \cite{gada}. We varied
the degree of smoothness by adjusting the \code{T} parameter of the
\code{BackwardElimination} function, and for the \model{gada.default}
model, we did not use the \code{BackwardElimination} function.

We used version 1.3 of the \package{flsa} package from CRAN to
calculate the Fused Lasso Signal Approximator
\cite{fused-lasso-path}. The FLSA solves the following optimization
problem for each chromosome:
\begin{equation}
  \label{eq:flsa}
  \hat x^\lambda = 
\argmin_{\mu\in\RR^m} 
\frac 1 2 \sum_{i=1}^d (x_i-\mu_i)^2
+\lambda_1\sum_{i=1}^d|\mu_i|
+\lambda_2\sum_{i=1}^{d-1}|\mu_i-\mu_{i+1}|.
\end{equation}
We define a grid of values $\lambda\in\{10^{-5},\dots,10^{12}\}$, take
$\lambda_1=0$, and consider the following parameterizations for
$\lambda_2$:
\begin{itemize}
\item \model{flsa}: $\lambda_2=\lambda$.
\item \model{flsa.norm}: $\lambda_2=\lambda m \times 10^6/l$ where
  $m$ is the number of points and $l$ is the length of the chromosome
  in base pairs.
\end{itemize}


We used version 1.18.0 of the \package{DNAcopy} package from
Bioconductor to fit a circular binary segmentation model
\cite{dnacopy}. We varied the degree of smoothness by adjusting the
\code{undo.SD}, \code{undo.prune}, and \code{alpha} parameters of the
\code{segment} function.  However, the dnacopy.prune algorithm was too
slow ($>24$ hours) for some of the profiles with many data points, so
these profiles were excluded from the analysis of dnacopy.prune.

We used version 0.2.1 of the \package{cghFLasso} package from CRAN,
which implements a default fused lasso method \cite{cghFLasso}, but
does not provide any smoothness parameters for breakpoint detection.

We used version 2.0.0 of the \package{GLAD} package from Bioconductor
to fit the GLAD adaptive weights smoothing model \cite{glad}. We
varied the degree of smoothness by adjusting the \code{lambdabreak}
and \code{MinBkpWeight} parameters of the \code{daglad} function. For
the \model{glad.haarseg} model, we used the
\code{smoothfunc="haarseg"} option and varied the \texttt{breaksFdrQ}
parameter to fit a wavelet smoothing model \cite{haarseg}.

To fit a Gaussian maximum-likelihood piecewise constant smoothing
model \cite{statistical-approach}, we used pruned dynamic programming
as implemented in version 0.1 of the \package{cghseg} package from
R-Forge \citep{pruned-dp}. For the default \model{cghseg.mBIC} model,
we used the modified Bayesian information criterion 
\cite{mBIC}, which has no smoothness parameter, and is implemented in
the \code{uniseg} function of the \package{cghseg} package.  For the
\model{cghseg.k} model, we used the \code{segmeanCO} function with
\code{kmax=20} to obtain the maximum-likelihood piecewise constant
smoothing model $\mu^k\in\RR^m$ for $k\in\{1, \dots, 20\}$ segments.
Lavielle suggested penalizing $k$ breakpoints in a signal sampled at
$m$ points using $\lambda k$, and varying $\lambda$ as a tuning
parameter \cite{lavielle2005}. We implemented this model selection
criterion as the \model{cghseg.k} model, for which we define the
optimal number of segments
\begin{equation}
  \label{eq:cghseg.k}
  k^*(\lambda) = 
\argmin_{k\in\{1,\dots,20\}}
\lambda k+  \frac 1 m \sum_{i=1}^m (x_i-\mu_i^k)^2,
\end{equation}
and the optimal smoothing $
  \hat x^\lambda = \mu^{k^*(\lambda)}.
$

We used the \code{cpt.mean} function in version 1.0.4 of the
\package{changepoint} package from CRAN to fit a penalized maximum
likelihood model using a Pruned Exact Linear Time (PELT) algorithm
\cite{pelt}. PELT defines $\mu^k$ in the same way as cghseg, but
defines the optimal number of segments as
\begin{equation}
  \label{eq:pelt}
  k^*(\beta) = \argmin_{k\in\{1,\dots,m\}}
 \beta (k-1) + \sum_{i=1}^m (x_i-\mu_i^k)^2.
\end{equation}
For the \model{pelt.default} model, we used the default settings which
specify \code{penalty="SIC"} for the Schwarz or Bayesian Information
Criterion, meaning $\beta = \log d$. For the \model{pelt.n} model, we
specified \code{penalty="Manual"} which means that the \code{value}
parameter is used as $\beta$, and the \code{cpt.mean} function returns
$\mu^{k^*(\beta)}$. We defined the same grid of $\lambda$ values that
we used for \model{cghseg.k}, and let $\beta = \lambda d$. Note that
this model is mathematically equivalent to \model{cghseg.k}, but shows
small differences in the results, since there are rounding errors when
specifying the penalty \code{cpt.mean(value=sprintf("n*\%f",lambda))}
for \model{pelt.n}.








    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Authors contributions}
OD, GS, IJ-L, and JC designed and executed the experiments. TDH, VB,
GS, and IJ-L created the annotation databases. TDH, FB and J-PV
conceived the statistical model selection framework. TDH wrote the
code and manuscript. All authors read and approved the final
manuscript.

    

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgements}
\ifthenelse{\boolean{publ}}{\small}{} Thanks to Edouard Pauwels for
helpful comments to simplify the mathematics on an early draft of the
paper. Thanks to the two reviewers for encouraging us to show multiple
annotation data sets and a high-density profile.

This work was supported by Digiteo [DIGITEO-BIOVIZ-2009-25D to
T.D.H.]; the European Research Council [SIERRA-ERC-239993 to F.B; SMAC-ERC-280032 to J-P.V.];
the French National Research Agency [ANR-09-BLAN-0051-04 to J-P.V.]; the Annenberg Foundation [to G.S.];
the French Programme Hospitalier de Recherche Clinique [PHRC IC2007-09
to G.S.]; the French National Cancer Institute [INCA-2007-1-RT-4-IC to
G.S.]; and the French Anti-Cancer League.




 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%              
%%  Bmc_article.bst  will be used to                       %%
%%  create a .BBL file for submission, which includes      %%
%%  XML structured for BMC.                                %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %% 
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %% 
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


{\ifthenelse{\boolean{publ}}{\footnotesize}{\small}
 \bibliographystyle{bmc_article}  % Style BST file
  \bibliography{refs}
 }     % Bibliography file (usually '*.bib' ) 

%%%%%%%%%%%

\ifthenelse{\boolean{publ}}{\end{multicols}}{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files

\newpage

\section*{Figures}

\includegraphics[width=\textwidth]{figure-1-smoothing}

\subsection*{Figure 1 - Breakpoint annotations quantify the accuracy
  of a smoothing model}
Model agreement to annotated regions can be measured by examining the
positions of predicted breakpoints $\hat b_{i}^\lambda$ (vertical
black lines) observed in the smoothing model $\hat y^\lambda$ (green
lines).  Black circles show logratio measurements $y$ plotted against
position $p$ for a single profile $i=375$. Chromosomes are shown in
panels from left to right, and different values of the smoothing
parameter $\lambda$ in the \model{flsa} model are shown in panels from
top to bottom. Models with too many breakpoints ($\lambda=0.5$) and
too few breakpoints ($\lambda=10$) are suboptimal, so we select an
intermediate model ($\lambda=7.5$) that maximizes agreement with the
annotations, thus detecting a new breakpoint on chromosome 7 which was
not annotated.

\newpage

\includegraphics[width=\textwidth]{figure-2-learning-curves}

\subsection*{Figure 2 - Selecting the model that minimizes the
  breakpoint annotation error}
Training error functions for global and local models plotted against
smoothing parameter $\lambda$. The original annotation data set was used to
calculate the annotation error. In the top row panels, we plot
$E^{\text{global}}(\lambda)$ from Equation~\ref{eq:globalerror}, and
in the other rows, we plot $E^{\text{local}}_i(\lambda)$ from
Equation~\ref{eq:localerror}. Each column of plots shows the error of
a particular algorithm, and the minimum chosen using the global
training procedure is shown using a vertical grey line. Note that the
local model training error can be reduced by moving from the globally
optimal smoothing parameter $\hat \lambda$ to a local value $\hat
\lambda_i$, as in profile $i=375$ for dnacopy.sd. For the local models
trained on single profiles, many smoothing parameters attain the
minimum. So we use the protocol described in the ``\sectionpick''
section to select the best value, shown as a black dot.

\newpage

\includegraphics[width=\textwidth]{figure-ann-sets}

\subsection*{Figure 3 - Two annotation data sets of the same profile}

The annotated chromosomes of profile $i=8$, with the two sets of
annotations shown in the two rows of plots. For the original
annotations, the same 6 regions on each profile were systematically
labeled as either 0breakpoints or $>$0breakpoints. For the detailed
annotations, a GUI was used to draw regions anywhere on the profile,
and 1breakpoint annotations were also used.

\newpage

\includegraphics[width=\textwidth]{figure-roc2}

%\includegraphics[width=\textwidth]{figure-3-roc.pdf}

\subsection*{Figure 4 - ROC curves compare breakpoint detection of
  global models on two annotation data sets}
ROC curves for the training error with respect to the breakpoint
annotation data are shown as colored lines. The curves are shown in 3
panels zoomed to the upper left region of ROC space to avoid visual
clutter. Each curve is traced by plotting the error of a model as the
degree of smoothness is varied, and an empty black circle shows the
global model chosen by minimizing the error with respect to all
annotations. Algorithms with no tuning parameters are shown as black
dots. Note that some ROC curves appear incomplete since some
segmentation algorithms are not flexible enough for the task of
breakpoint detection, even though we ran each algorithm on a very
large range of smoothness parameter values.

\newpage
\includegraphics[width=\textwidth]{figure-global-local}

\subsection*{Figure 5 - Global models are better breakpoint detectors
  than local models}

Leave-one-out cross-validation over the annotions on each profile was
used to compare breakpoint detection error of global and local
models. The two panels show the two annotation data sets, and each row
shows the performance of one of the models described in the Algorithms
section. After selecting the smoothness parameter $\lambda$ by
minimizing either the global or local annotation error, we plot the
mean and standard deviation over 10 test sets. Each default model does
not have a smoothness parameter, and shows equivalent local and global
model error. Squares show the same colors as the other figures and
tables, and are absent for default models that have no smoothness
parameters.

\newpage

\includegraphics[width=\textwidth]{figure-6-learned}

\subsection*{Figure 6 - Learned global models are better breakpoint
  detectors than default models}
Chromosomes with detailed annotations for neuroblastoma copy number
profile $i=207$ are plotted in panels from left to right, and each row
shows one of 5 segmentation models drawn in green. Breakpoints
detected by the models are drawn as vertical black lines. The bottom 3
rows show default models trained using no annotations, and the top 2
rows show global models trained using annotations of other
profiles. The equivalent cghseg.k and pelt.n global models show only 2
annotation errors compared to 3 and 10 in the corresponding default
models. Furthermore, the equivalent cghseg.k and pelt.n global models
perform better than dnacopy.sd, which shows 4 annotation errors.

\newpage



\includegraphics[width=\textwidth]{figure-test-error-train-profiles}

%\includegraphics[width=0.8\textwidth]{figure-5-kinetics}

\subsection*{Figure 7 - Test accuracy increases to a model-specific
  limit as number of training profiles increases}
Cross-validation was used to estimate the generalization ability of
the global models with different sized training sets. Panels from left
to right specify the annotations that were used to train, and panels
from top to bottom specify the annotations that were used to
test. Note that the results change very little even when training on
one data set and testing on another. For each training set size $t$,
the profiles were partitioned into training sets of approximately size
$t$, then were evaluated using the annotations from all the other
profiles. Results on these data indicate increasing accuracy (lines)
and decreasing standard deviation (shaded bands) as the training set
increases. The accuracy of each model quickly attains its maximum,
after only about $t=10$ profiles. In Table~4, we show the results for
all algorithms when trained on $t=10$ profiles from the detailed data
set, and tested on the other profiles in the detailed data set
(vertical black line).

\newpage

\includegraphics[width=\textwidth]{figure-snp6-both}

\subsection*{Figure 8 - Zooming allows creation of small annotated
  regions on high-density profiles}

\textbf{Top}: three regions of chromosome on from a high-density
Affymetrix SNP6 array is shown along with some breakpoint
annotations. There are almost 2 million probes on the array,
\input{snp6-chr2-probes}\unskip\ probes on chromosome 2, and
\input{snp6-chr2-shown}\unskip\ probes shown in these 3
windows. \textbf{Bottom}: zoom to show detail of windows 2 and
3. Annotations can be used to ensure that a smoothing model accurately
recovers breakpoints around small $\approx 10$kb alterations such as
those shown in windows 2 and 3.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage

\subsection*{Table 1 - Counts of annotations in two annotation data
  sets of the same copy number profiles}

\begin{center}
  \input{table-data-sets}
\end{center}

For the two annotation data sets (columns), we show the annotation
protocol, counts of annotated profiles and chromosomes, and counts of
annotations.

%% Use of \listoftables is discouraged.
%%
% \section*{Tables}
% \subsection*{Table 1 - Counts of annotations in the neuroblastoma data
%   set conditional on region}
% \input{table-annotation-chrom-counts}

% Counts of normal and breakpoint annotations in the neuroblastoma data
% set, conditional on region. Min and max limits of each region are
% shown in mega base pairs, in reference to the Hg19 Human genome
% assembly.

% \subsection*{Table 2 - Counts of profiles in the neuroblastoma data set
%   conditional on annotations}
% \input{table-annotation-profile-counts}

% Counts of profiles in the neuroblastoma data set, conditional on
% number of annotations.  Note that most profiles have more normal
% regions than breakpoint regions.  For example, 335 profiles have all 6
% regions annotated as normal.

%\subsection*{Table 3 - Estimated breakpoint detection error on
%  un-annotated regions}
%\input{table-generalization-error-global-models}

\subsection*{Table 2 - Training error of local models}

\begin{center}
  \input{table-training-error}
\end{center}

For both the original and the detailed annotation databases, and for
all 17 algorithms, we report the mean training error of the local
model in percent. For each profile and algorithm, the smoothing
parameter with minimal breakpoint annotation error was selected, and
we report the mean training error across all profiles.

\newpage

\subsection*{Table 3 - Ranking global models by breakpoint detection
  error on un-annotated profiles}
%\input{table-error-on-unseen-profiles}
\begin{center}
  \input{table-unseen-detailed}
\end{center}

The $n/t$-fold cross-validation protocol was used to estimate error,
false positive (fp), and false negative (fn) rates in the detailed
annotations of the $n=575$ profiles. Squares show the same colors as
in the figures, and are absent for default models that have no
smoothness parameters. The smoothness parameter was chosen using
annotations from approximately $t=10$ profiles, and mean and standard
deviation (sd) of test error over $\lfloor n/t\rfloor =57$ folds are
shown as percents. The Timings column shows the median number of
seconds to fit the sequence of smoothing models for a single profile.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\end{bmcformat}
\end{document}







